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1 Introduction 

The most widely known and used result in asymptotics is probably Stirling 's 
formula, 

n\ ~ V27rnn"e"", n ^ oo (1) 

named after James Stirling (May 1692 - 5 Dec. 1770). The formula pro- 
vides an extremely accurate approximation of the factorial numbers for large 
values of n. The asymptotic formula for which Stirling is best known, 
appeared as Example 2 to Proposition 28 in his most important work, Metho- 
dus Dijferentialis, published in 1730 [Hj. In it, he asserted that log(n!) is 
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approximated by "three or four terms" of the series 

n+01og(^n + _a(^n + +ilog(27r) (2) 

a 7a 
+ 



24(n+i) 2880 (n+i)' 

where log means the base-10 logarithm and a = [ln(10)]~^ [HI] . 

In 1730 Stirling wrote to Abraham De Moivre (26 May 1667 - 27 Nov. 
1754) pointing out some errors that he had made in a table of logarithms of 
factorials in the book and also telling him about (j21). After seeing Stirling's 
results, De Moivre derived the formula 

1„ |(„ _ . („ _ ,„(„, (2.) + (3) 



^2A; {2k 



which he published in his Miscellaneis Analyticis Supplementum a few months 
later. Equation (jHl) is called Stirling's series and the numbers Bk are called 
the Bernoulli numbers, and are defined by 



Bo=l, E(^)^)i?.=0, k>l. (4) 

Clearly Stirling and De Moivre regularly corresponded around this time, for 
in September 1730 Stirling relates the new results of De Moivre in a letter 
to Gabriel Cramer. 

In 1729 Leonhard Euler (15 April 1707 - 18 Sept. 1783) proposed a 
generalization of the factorial function from natural numbers to positive real 
numbers ^I]. It is called the gamma function, T{z), which he defined as 

r(2;) = hm . . r , \ ^ (5) 

n^ooz[Z + 1) ■ ■ ■ [Z + n) 

and it is related to the factorial numbers by 

r(n + l) = n!, n = 0,1,2,.... 

From Euler's definition (0), we immediately obtain the fundamental relation 

T{z + 1) = zT{z) (6) 

and the value r(l) = 1. In fact, the gamma function is completely charac- 
terized by the Bohr-Mollerup theorem |HI: 
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Theorem 1 The gamma function is the only function T : (0, oo) (0, oo) 
which satisfies 

1. r(i) = 1 

2. r(a;+ 1) = xT{x) 

3. In [r(x)] is convex 
for all X G (0, oo) . 

Proof. See g] and ■ 

Another complex-analytic characterization is due to Wielandt |4Hj : 

Theorem 2 The gamma function is the only holomorphic function in the 
right half plane A satisfying 

1. r(l) = 1 

2. V{z + 1) = zV{z) for allzeA 

3. r (z) is bounded in the strip 1 < Re(2;) < 2 

Proof. See [35j. ■ 

In terms of T{z), we can re- write as 

\n[T{z)]r^r{z), z^oo (7) 

with 

Viz) = zln{z) - z-]-ln{z) + ^ln{27v) (8) 
and in the form 

\n[T{z)]^V{z)+RN{z), z^oo (9) 
where Ro{z) = and 

^ B 

Estimations of the remainder In [r(2;)] —V{z) — R]\f{z) were computed in 



2 Previous results 



Over the years, there have been many different approaches to the derivation 
of ((Tj) and Q, including: 

1. Aissen ^ studied the sequence Vn = " ■ Using his lemma 
Lemma 3 // 

^ = 1 + ^ + O (n-2) 
Vn n ^ ' 

and yn ^ for all n, then i/n ~ Cn", n —* oo for some non-zero 
constant C. 

he showed that n\ ~ C ^/nvH^ e'"^ . 

2. Bender & Orszag Bleistein & Handelsman [Hj, Diaconis & Freedman 
[ini, Dingle [Hj, Olver [SS| and Wong [IS] applied Laplace's method to 
the Euler integral of the second kind 

oo 

T{z)= j e-\-'dt, Re{z) > 0. (11) 



3. Bender & Orszag j5j and Temme PP] used HankeVs contour integral [T] 

(0+) 

= — [ t^'e'dt 
Tiz) 27ri J 

— oo 

and the method of steepest descent. 

4. Bleistein & Handelsman [Hj, Lebedev j2S], Sasvari and Temme [10] 

used Binet's first formula 

oo 

In [r(^)] = P(^) + |i Q - 1 + e-*^dt, Re(^) > 



and 

t\2 t e*-i; ^(2A;)! " 
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5. Blyth & Pathak T] and Khan [221 used probabilistic arguments, ap- 
plying the Central Limit Theorem and the limit theorem for moment 
generating functions to Gamma and Poison random variables. 



6. Coleman P defined 



n -\ — ) ln(n) — n + 1 — In (n!) , 



and showed that Cn —>■ 1 — | In (27r) as n — > oo. A similar result was 
proved by Aissen |2], using the concavity of ln(a;). 



7. Dingle ^| used Weierstrass ' infinite product 

oo 



^ oo 

1 ^^T-r r/. Z 



(where 7 is Euler's constant) and Mellin transforms. 
8. Feller [12], [IB] proved the identity 



1 fl\ '"'^ 

In (n!) - - ln(n) = I{n) - / ( - j + ^ (a, - 6,) + a„, (12) 

^ ^ fc=i 



where 

I{n)= [\n{t)dt, ak= [ \n(-]dt, bk = [ \n ( y ] dt 



t \k 



and showed that 



There is a mistake in his Equation (2.4), where he states that 

1 /1\ """^ 

In (n!) - - ln(n) + I{n) " ^ ( 2 ) = ^(^k - h) + an, 

instead of (|T5|) . 
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9. Hayman ^\ used the exponential generating function 



k=0 

and his method for admissible functions. 
10. Hummel established the inequalities 

^ <r„ + iln(27r) < 1, n = 2,3, 



where 



r„ = In 



n\e 



V2 



Tcnn' 



Impens [21], showed that for x > 

R2n{x) < In [r(x)] - V{x) < R2m+l{x), 71,171 >0, 

where Rn{x) was defined in (fTUj). Maria j26j showed that 

3 



12n + ■ 

2(2n + l)_ 

Mermin j2Hl proved the identity 



< r„, n=l,2,.... 



k=n 



which he used to show that r„ ~ Rsin). Michel f2^, proved the in- 
equality 



e"" - 1 



1 



1 



< 



1 



1 



360n3 108n4 ' 



12n 288n2 
Nanjundiah showed that 

R2{n) < < Ri{n) = 1,2, . . . . 
Robbins [36j, established the double inequality 

< r„, < n = 1,2, . . 



n = 3,4.... 



12n+ 1 



12n 



11. Marsaglia & Marsaglia [27] derived from (fTTj) the asymptotic expansion 

„!^„-e-»|.,(l)*r(|)., 

where the generating function G{z) = ''^^b^z'' is defined by 

G{z) exp [1 - G{z)] = exp (^-^^^) . ^'(0) = 1- 

12. Namias introduced the function F{n) = with V{n) defined in 
(jHI). From Legendre's duplication formula 

r(2n) = — ^r(n)r(n + -j, (13) 



he derived a functional equation for F{n) 



F{n)F{n-l) \ 2n 

from which he obtained ©• He also considered the triplication case, 
using Gauss' multiplication formula 

Y{mz) = {^p2^ ""m"'"-^ JJr ( 2 + — j , m = 2,3,... 

with m = 3. His results where extended by Deeba & Rodriguez in ^2)- 

13. Olver used Euler's definition (0) 

In [r(z)] = - In (2) + lim z In (n) + [In (A;) - In (2 + fc)] 

n— >oo ^— ^ 

fc=l 

and the Euler-Maclaurin formula. A similar analysis was done by 
Knopp 1211 and Wilf -44. 

14. Patin used ()11|) and the Lebesgue Dominated Convergence Theo- 
rem. 
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15. Whittaker & Watson |12| used Binet's second formula 



i Re(z)>0 



and 

arctania;) = xy — 

k=0 



\x\ < 1. 



Thus, there have been a huge variety of approaches to Stirhng's result, 
ranging from elementary to heavy-machinery methods. In an effort to join 
such illustrious company, we present still another direction for deriving ©. 

Our starting point shall be the difference equation (jHl). A parallel ap- 
proach was considered in [Sj. For a different analysis of © using the method 
of controlling factors, see jS]. Extensions and other applications of the 
method used can be found in jTU], [T3] and [TB] . 

3 Asymptotic analysis 
3.1 Stirling's formula 

We begin with a derivation of ((Zj), to better illustrate how the method works. 
We assume that 

\n[T{z)]r^ f{z)+g{z), z ^ oo (14) 

with 

9 = oif), z-^oo. (15) 

Using (fT^ in 0, we have 

f{z + 1) - f{z) + g{z + 1) - g{z) ~ \n{z). (16) 
Expanding f{z + 1) and g{z + 1) in a Taylor series, we obtain 

f\z) + ]^f"{z)+g'{z)r^\n{z). (17) 

From ()15|) and p7|) we get the system 

nz)=Hz), ]^f{z)+g'{z) = Q 
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and thus, 

f{z) = z\n{z)-z, g{z) = -^\n{z) + C. (18) 
To find the constant C in |T8|) . we replace T{z) ~ e-^^^^^^^^^ in (|T3|l and obtain 

222-1 



J{2z)+a{2z) ^ ^!lJ.g/(2)+<?Wg/(^+i)+3(^+|)_ 



or 



e^~^ + — j ~ V27r, ^ ^ oo, 

from which we conclude that C = In (v^27r) . 
Hence, we have shown that 

1 1 

In [r(2;)] ~ z \Yi.{z) ~ ~ 2 ^^('^) ~'~ 2 ^'^'^'^ ' ~^ 
3.2 Stirling's series 

To extend the result of the previous section, we now assume that 

N 

ln[r(;2)] z ^ 00 (19) 

fc=0 

with 

/fc+i = o(/fc), ^^oo, A; = 0,l,...,iV-l. (20) 
Using (fT^ in © we have 

^/fc(z + l)-/fe(z)~ln(z), z^^. (21) 

fc=0 

Replacing the Taylor series 

+ 1) = E^IjA(-) 

i>0"'' 



in (PT|). we have 



^ 1 rf^- 



EE7y^/^(^)-i^(^)' (22) 



fc=0_ 
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From (j2(Jj) . we obtain the system of ODEs 



and 

which imply 

foiz) = z\n{z) - z (23) 

and 

k 

^(-) = ^^1' (24) 

where we have omitted (for the time being) any constant of integration. To 
find the functions fk{z), we set 

fkiz)=ak^foiz) (25) 

in and get ao = 1 and 

d^ ^ 1 d^ d^~^ 

which gives 

k ^ 

ao = 1, ak = -J^^— k>l, 

or 



ao = l, Et^-T— = 0' ^^1- (26) 



3=0 

Multiplying both sides by (A; + 1)!, we can write (pUj) as 



^.^J/c + l-j)!j 



10 



or 

ao = l, X^(^^^)j!ai = 0, k>l. (27) 

Comparing (jD) and (j77j) we conclude that 

afc = ^, A; > 0. (28) 
Thus, from (|2S1), (|2S1) and (EH) we have 

from which we obtain 

A(z) = -iln(z) (29) 

and 

Since i?2fc+i = for all /c > 1, we need to consider even values of k only, 

So far, we haven't included any constant of integration in our calculations. 
We could add a constant to one of the functions fk{z), let's say to fi{z), and 
proceed as in Section I^TTl to find it. Doing this, we would obtain from ()2H|1 . 
(j^ and dnni) that 

N 
k=0 

where V{z), Rn{z) were defined in (jSJ and (fTUj) respectively. 

Another possibility, would be to assume no previous knowledge of T{z), 
except for the difference equation r{z + 1) = zr{z) and the value at 1, 
r(l) = 1. In doing so, (fTIH) would imply the initial conditions /fc(l) = 0, for 
all A; > 0. Hence, we would have 

1 

fo{z) = z\n{z) - z + 1, fi{z) = --ln{z) 

~ 2k{2k + l)z^>'-^ 2k{2k+l)' - 
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and therefore 

In [T{z)] ~ z \n{z) - ^ -\ ^^{A + C^tv + Rn{z), (31) 

with 

Computing the first few Cat, we would get 

Ci ~ .91667, C2 ^ .91944, C3 ~ .91865, C4 ~ .91925, C5 ~ .91840, 
~ .92032, C-j - .91391, Cg ~ .94346, Cg ~ .76382, Cio ~ 2.1562, 

and increasingly greater numbers (in absolute value). We would conclude 
that, before the sum starts diverging, the Cat seem to be approaching a 
value close to .918. Given our previous discussion of Stirling's formula, it is 
not surprising to find that \ In (27r) ~ .91894. Of course, geniuses like Euler or 
Gauss would reach the conclusion that the optimal constant equals | In (27r) 
without knowing anything about Stirling's work! 

4 Conclusion 

We have presented the history and previous approaches to the proof of Stir- 
ling's series dH}. We have used a different procedure, based on the asymptotic 
analysis of the difference equation (jHl). The method reproduces dH)) very easily 
and can be extended to use in more complicated difference equations. 
Bender and Orszag observed in ^5, Page 227] that 

without further information the constant \ In (27r) cannot be 
determined. The difference equation that we have solved is linear 
and homogeneous, so any arbitrary multiple of a solution is still 
a solution. 

While agreeing with them completely, we have shown that by imposing 
the additional condition r(l) = 1 one can find an approximation to the value 
of I In (27r) , without any other assumptions. Thus, local behavior at 2 = 1 
and asymptotic behavior as z — > cxd can be combined fruitfully. 

We sincerely hope that more and more professors and their students will 
discover the beauty contained in the (very!) Special Functions, among which 
r(2;) is, without doubt, a prima donna. 
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